Density minimum and liquid-liquid phase transition 
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We present a high-resolution computer simulation study of the equation of state of ST2 water, 
evaluating the liquid-state properties at 2718 state points, and precisely locating the liquid-liquid 
critical point (LLCP) occurring in this model. We are thereby able to reveal the interconnected set 
of density anomalies, spinodal instabilities and response function extrema that occur in the vicinity 
of a LLCP for the case of a realistic, off-lattice model of a liquid with local tetrahedral order. In 
particular, we unambiguously identify a density minimum in the liquid state, define its relationship 
to other anomalies, and show that it arises due to the approach of the liquid structure to a defect-free 
random tetrahedral network of hydrogen bonds. 

PACS numbers: 61.25.-f,64.30.+t,64.70.-p 



Water is not only the most abundant liquid on earth, 
but also a prototype of many network forming materials. 
As for water, important substances such as silicon and 
silica exhibit in their liquid phase a disordered network 
structure that arises due to highly directional tetrahe- 
dral bonding interactions. It has long been appreciated 
that the development of an open tetrahedral network on 
cooling is related to the occurrence of the density maxi- 
mum observed in these liquids 0. In water the density 
maximum occurs at 277 K at ambient pressure P. For 
temperatures T above the density maximum, the isobaric 
expansivity ap is positive, as for the vast majority of liq- 
uids. As T decreases through the density maximum, ap 
decreases to zero and then becomes negative; in this re- 
gion the liquid expands on cooling. For all liquids with 
density maxima, an open question concerns what hap- 
pens as T decreases further: Does ap remain negative, 
or eventually return to "normal" positive values? In the 
latter case, the T at which ap again becomes positive 
corresponds to a density minimum. 

A density maximum, although well-known in the case 
of water, occurs in only a few liquids. Experimental ob- 
servations of liquid-state density minima are almost non- 
existent; we have found only one report in the literature, 
a study of a mixture of Ge and Se 0- Yet the possibil- 
ity of a density minimum occurring in deeply supercooled 
water and similar systems has been discussed in the liter- 
ature (see e.g. Ref. Q), and is supported by several theo- 
retical and simulation studies. Density minima have been 
identified in lattice models of network- forming fluids , 
and recent computer simulations of water also 

present data in which a liquid-state density minimum, 
or a trend toward one, is seen or can be inferred. Also 
notable is the recent experimental report of a density 
minimum in vitreous silica, although non-equilibrium ef- 
fects complicate the interpretation of the results in terms 
of equilibrium liquid-state behavior 

At the same time, a growing body of evidence sug- 
gests that a liquid- liquid critical point (LLCP) may also 
be a characteristic feature of liquids with local tetrahe- 



dral order [1 H EH, El . Below such a LLCP, two dis- 
tinct liquid phases, the so-called "low density" and "high 
density" liquids (respectively LDL and HDL) are sepa- 
rated by a line of first order phase transitions. Formal 
thermodynamic analysis has revealed general relation- 
ships that must be satisfied in a liquid displaying both 
density anomalies and the anomalies associated with a 
LLCP [1 [H Q E [S E3- Some studies, such as 
in Ref. |8j, associate multiple density anomalies with a 
polyamorphic transition related to a LLCP. However, an 
unambiguous realization of these relationships for a re- 
alistic, off-lattice, equilibrium molecular liquid has not 
been reported, either experimentally or in computer sim- 
ulations. In particular, the implications of a liquid-state 
density minimum merit focused study, both in terms of 
its relationship to the LLCP, as well understanding its 
origins in terms of the low T thermodynamic and struc- 
tural behavior of the liquid. For example, if a density 
maximum is followed at lower T by a density minimum, 
it means that many, if not all of the anomalies for which 
water in particular is famous, ultimately disappear under 
sufficient supercooling. The possibility of such a return 
to "normal" behavior has important implications for un- 
derstanding the properties of bulk supercooled water, the 
amorphous ices, as well as composite systems in which 
the tetrahedral network is influenced by surfaces or so- 
lutes H0. 

In this Letter, we present a high- resolution computer 
simulation study of the equation of state (EOS) of ST2 
water E3 , a molecular dynamics model known to exhibit 
a prominent density maximum and also a clear example 
of a LLCP 9]. As shown below, we find in addition 
that this model liquid displays a density minimum. We 
are thereby able to reveal the thermodynamic context in 
which this density minimum occurs, identify the struc- 
tural origins of the anomaly, and clarify its relationship 
to the LLCP. 

Our results are based on molecular dynamics simula- 
tions of liquid systems of N = 1728 molecules. All simu- 
lations are carried out at constant N and volume V, and 
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T is maintained at the desired value using Berendsen's 
method |19| with a time constant for the relaxation of T 
fluctuations of 2 ps. Molecular interactions are cut off 
at a distance of 0.78 nm, and the reaction field method 
is used to approximate the contribution of dipole-dipolc 
interactions at longer range. The time step used is f fs. 

Our simulation protocol is as follows: Randomized 
starting configurations are prepared at each density p 
to be studied and are used to initiate runs at 400 K. An 
equilibration and production phase are carried out, after 
which the target T is reduced by 5 K, and a new equili- 
bration/production cycle is initiated using the last con- 
figuration of the previous run. In this way, the properties 
of successive state points along the specified isochore (i.e. 
at constant p) are evaluated with decreasing T. Both the 
equilibration and production phases at each T are carried 
out for the time required for the mean squared displace- 
ment (MSD) to reach 1 nm 2 , or for 100 ps, whichever 
is larger. Note that a MSD of 1 nm 2 corresponds to 
an average displacement of about 10 times the distance 
between nearest neighbour molecules in liquid water at 
p = 1.0 g/cm 3 . For our lowest T state points, our run- 
time criterion leads to simulations of on the order of 
100 ns. Using the above protocol, we evaluate the prop- 
erties of the liquid at 2718 state points: at densities from 
p = 0.70 to 1.50 g/cm 3 in steps of 0.01 g/cm 3 ; and for T 
from 400 K down to at least 255 K, in 5 K steps. Along 
some isochores, a T as low as 200 K is reached. The sum 
of the computational work reported here corresponds to 
approximately 30 cpu- years on a 3.0 GHz Intel Xeon pro- 
cessor. 

The resulting EOS data for P(p, T) is shown in 
Fig. nja). Early support for the possibility of a LLCP 
in water was based on the shape of the ST2 EOS re- 
ported in Ref. 9|. These features of the ST2 model are 
seen here in much greater detail. The LLCP is readily 
located, as well as the metastability limits (spinodals) of 
both the LDL and HDL phases that meet at the critical 
point. We also identify an upper bound on the location 
of the liquid-gas spinodal by locating the cavitation limit 
of the liquid at the extreme low p range at negative P. 

Fig.^a) shows most of the individual state points sim- 
ulated, as well as isochores obtained from the fit of an 
analytic EOS to the data set [2(j. Points on the locus of 
density maxima A max are coincident with the minima of 
these isochores. As p decreases, A max follows a retracing 
path in the T-P plane to p — 0.82 g/cm 3 , but notably, 
the isochores for p = 0.81 and 0.80 g/cm 3 do not ex- 
hibit minima. Rigorous thermodynamic arguments have 
shown that A max cannot end at finite T in isolation from 
other anomalies [l3L ITiL fl5| . It must either terminate 
on a spinodal locus, or A max itself must pass through an 
extremum in the T-P plane, at which point it becomes 
a line of density minima, A m j n . We find that the latter 
possibility occurs, as shown by both the shape of A max 
and the EOS isochores: low T maxima in the p = 0.82 
and 0.83 g/cm 3 isochores correspond to points on A m i n . 




-200 - 



200 225 



250 275 
T(K) 



300 325 350 



-70 
-80 
-90 
£ -100 
~-110 
-120 
-130 
-140 





i6- 




3*5 2 




- 












— 


□V 














■ oX X 


~ ° n 


00 — 








• 0.82 g/cm 3 , N=216 X 




■ 0.83 g/cm 3 , N=216 - 




O 0.82 g/cm 3 , N=1 728 ; 


: (b) 


□ 0.83 g/cm 3 , N=1 728 z 


: , , 9 l , , , l , , 


t I i i i I i i i I i i i I i - 



240 260 280 



300 320 
T(K) 



340 360 



To further confirm the existence of A n 



we conduct 



FIG. 1: (color online) (a) Phase diagram of the N = 1728 
ST2 liquid. State points simulated are shown as small dots. 
Isochores (thin lines through dots) obtained from our fitted 
EOS are shown from p = 0.80 g/cm 3 (bottom- most line) and 
greater in steps of 0.01 g/cm 3 . Also shown is the estimated 
location of the liquid-gas spinodal (diamonds) , the LDL spin- 
odal (up-triangles), the HDL spinodal (down-triangles), and 
the LLCP (circle). A max (thick line) transforms into A m i n 
(thin line) at low P. Squares locate points on A m i n obtained 
from (b). (b) Variation of P with T along the p = 0.82 and 
0.83 g/cm 3 isochores for both the N = 1728 and N = 216 
systems. For the N = 216 data, error bars are drawn at twice 
the standard deviation of the means of P as evaluated from 
each of the 40 runs conducted at each state point. 



a new set of simulations focusing on the p = 0.82 and 
0.83 g/cm 3 isochores, with the aim of reducing the statis- 
tical errors. These simulations employ a smaller system 
size, N = 216, but we average over an ensemble of 40 in- 
dependent runs for each state point |2l| . The variation of 
P with T along these two isochores [Fig. QJb)] confirms 
that both a density maximum and a density minimum 
occur. The locations of the density minima so identi- 
fied are shown as squares in Fig.^a), and are consistent 
with the predictions of the data set obtained using 1728 
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FIG. 2: (color online) Location of the LLCP (C), spinodals 
(triangles and diamonds), density extrema (A), and response 
function extrema (A and T) in (a) the T-P and (b) the p-T 
planes. The symbols used for the various spinodals are the 
same as in Fig.0 Thick lines indicate maxima, and thin lines 
locate minima. 
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FIG. 3: (color online) Properties of the liquid along the 
p = 0.83 g/cm 3 isochore as found from the N — 216 simula- 
tions, (a) Variation of U — 3RT with T for the liquid (circles), 
compared with U — 3RT for ice 1^ (line). R is the gas con- 
stant. The dashed line shows U - 3RT from the fitted EOS. 
U for ice 1^ is found from simulations of a proton-disordered 
configuration of iV = 432 ST2 molecules at p = 0.83 g/cm 3 . 
(b) The fraction of network defects fd as a function of T. fd 
is estimated by evaluating the equilibrium average fraction of 
molecules not having exactly four nearest neighbour O atoms 
within 0.35 nm, the distance of the first minimum in the O-O 
radial distribution function, (c) The skewness of the P(U) 
distributions plotted in Fig. 2] The skewness is evaluated as 
(1/n) X^iLiK^i ~~ U)/a] 3 , where n is the number of Ui values 
sampled, U is the average value of U, and a is the standard 
deviation of P(U). 



molecules. 

All along the locus A (the union of A max and A m ; n ) 
the condition (dp/dT)p — is satisfied. Formal thermo- 
dynamic analysis has shown that changes of sign of the 
slope of A in the T-P plane are associated with intersec- 
tions with certain response function extrema. Ref. fl6| 
shows that the point A in Fig. Ola), where A has infi- 
nite slope, is coincident with a point on a locus A along 
which (dKx/dT)p = 0, where Kt is the isothermal 
compressibility. By evaluating the appropriate deriva- 
tive of our fitted EOS, we confirm that this requirement 
is satisfied. We also identify a new thermodynamic rela- 
tion, derived by analogy to those in Refs. and 
that states that the point B in Fig. |2Ja,), where A has 
zero slope, is coincident with a point on a locus T along 



which (dCp / dP)x = 0, where Cp is the isobaric spe- 
cific heat. This requirement is also satisfied by our fitted 
EOS. Fig. HJb) shows the location of the features of (a) 
in the p-T plane. 

At the same time, the curves A and T are formally 
related to the LLCP. At a LLCP, denoted C, isotherms 
of Cp and isobars of Kt will both exhibit a divergence. 
For T > Tc, A and T radiate outward from C as lines 
of maxima that are the supercritical extension of the di- 
vergence occurring as T — > Tc- In the present data, 
these maxima transform into minima, and then intersect 
the line A, controlling the sign of its slope. While the 
existence of these response function extrema is not a suf- 
ficient condition for there being a LLCP at lower T, it 
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FIG. 4: (color online) Probability distributions of U along 
the p — 0.83 g/cm 3 isochore as found from the N = 216 
simulations. From left to right are shown the distributions 
for T = 255, 275, 295, 315, 335 and 355 K. R is the gas 
constant. 

is a necessary condition 0|. The interconnected pat- 
tern of density anomalies and response function extrema 
shown here therefore illustrates how the discovery of such 
features in experimental data can provide a rational pro- 
cedure for locating a LLCP at lower T. 

Examination of the thermodynamics and structure of 
the liquid in the vicinity of the density minimum reveals 
why it occurs. As shown in Fig. [31 both the potential en- 
ergy C7, and the fraction of network defects, fd, are at the 
lowest T behaving more and more as they would for a sys- 
tem approaching a random tetrahedral network (RTN) 
in which all nearest-neighbour intermolecular bonds are 
satisfied: the T dependence of U approaches an ice-like 
form, and fd approaches zero. The occurrence of inflec- 
tions in both these curves shows that the region of most 
rapid change has passed, and that the low T region in 
which the density minimum occurs corresponds to the 
final asymptotic consolidation of the liquid structure to- 
wards that of a perfect, defect-free RTN. It is plausible 
that cep should be positive (i.e. normal) for a defect-free 
RTN; e.g. cep > for crystalline ice I/j (an ordered tetra- 
hedral network) in this T range 22] . In this light, the 



density minimum simply reflects the fact that the pro- 
cess of RTN formation is coming to an end, and so the 
anomalous T dependence of thermodynamic properties 
returns to normal. 

Another indicator that the liquid is approaching the 
RTN state can be seen in the distribution of values of U, 
the potential energy, sampled in our simulations (Fig.0J. 
For T less than the inflection in U, the distribution P(U) 
narrows, and also loses its symmetrical shape and be- 
comes skewed, as quantified in Fig.^c). This narrowing 
and skewing suggests that there is a lower bound on the 
value of U, specifically U for the RTN, and that the sys- 
tem is entering the T regime where the existence of this 
lower bound influences the average properties |23j . 

The location and slope of the A m j n line in Fig. praises 
the possibility that a density minimum occurs in ST2 
water at positive P at sufficiently low T. This, and the 
generic nature of the mechanism underlying the density 
minimum found here, demonstrates the potential for the 
occurrence of density minima in the low T, P > be- 
haviour of any tetrahedral liquid having a density max- 
imum. More broadly, our results add to and illuminate 
the complex behaviour of tetrahedral liquids: The emer- 
gence of a RTN structure in these liquids can be under- 
stood as the source of their anomalies. Yet the eventual 
outcome of this RTN-forming process is to bring the liq- 
uid closer and closer to a set of near-perfect RTN states 
that lie at the bottom of the liquid's potential energy 
landscape jSJ 0, ■ Our simulations show that such a 
"liquid at the bottom of the landscape" represents a dis- 
tinct regime of behaviour: one that exhibits a number of 
exotic properties not observed at higher T (e.g. a crystal- 
like heat capacity) , and also one in which thermodynamic 
anomalies, so often associated with networking-forming 
liquids, disappear. 
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